%% Forecasts 
% Generate solution matrix
[GG,RR,CC,eu,SDX,ZZ,~,~,~,~,stanames,shockStru]=feval(model.handle,model.param,solveopt,addsol);
model.GG = GG;
model.RR = RR;
model.CC = CC;
model.eu = eu;
model.SDX = SDX;
model.ZZ = ZZ; 

%% Forecasts for observables 

% Extract forecasts for states and plot them 
nForc = (forecast.endForc-sampleStru.sampleVec(end))/0.25;

forecast.states={'Consumption', 'Hours',    'utilization (u)',    'k physical (k)',    'Investment','GDP'}; 
forecast.names ={'C detrended', 'hours',    'utilization'    ,    'k physical',  'I detrended'     ,'GDP detrended'}; 
forecast.location=cellposition(forecast.states,names.states); 
forecast.plot.rows=10; %number of rows per sheet
forecast.plot.cols=2; %number of columns per sheet

forecast.forcVec = (sampleStru.sampleVec(end):.25:forecast.endForc)';
forecast.forcsampcell=sample2date(forecast.forcVec);

yoffset=0.5; %used to make plots prettier 

[stateForc,obsForc]= forecast_loop(GG(:,:,end),RR(:,:,end),ZZ(:,:,end),CC(:,end),KFStru.smoothSt(end,:),[],nForc); 
tempTop=KFStru.smoothSt(end,posStru.obsInStates)+(ZZ(:,:,end)*CC(:,end))'; 
obsForc=[tempTop;obsForc]; 

%% Get States to Forecast
forecast.obs = { 'yobs (dy)' 'cobs (dc)' 'iobs (di)' 'Hours' 'FFR' 'JCXFEBM_LD100'};
forecast.scale = [4 4 4 1 4 4]';
forLocx4    = ismember(names.observables,{'yobs (dy)' 'cobs (dc)' 'iobs (di)' 'FFR' 'JCXFEBM_LD100'});
forLoc      = ismember(names.observables,{'Hours'});
tempobsForc = obsForc;
tempobsForc(:,forLocx4) = tempobsForc(:,forLocx4)*4;
allLoc      = forLocx4 | forLoc;

%% Forecast Obs
tempobsForc_Small = tempobsForc(:,allLoc');
figure ;
numForc=numel(forecast.obs);
forecast.obsSS = forecast.scale .* model.CC(cellposition(forecast.obs,names.states));

for j=1:numForc;
    subplot(ceil(numForc/forecast.plot.cols),forecast.plot.cols,j)
    plot(forecast.forcVec,tempobsForc_Small(:,j),'LineWidth',2)
    axis tight;
    
    fmin = min(min(tempobsForc_Small(:,j)),forecast.obsSS(j));
    fmax = max(max(tempobsForc_Small(:,j)),forecast.obsSS(j));
    ylim([fmin-yoffset fmax+yoffset]);
    
    hand = hline(forecast.obsSS(j));
    set(hand,'Color',[0.55 0.55 0.55],'LineWidth',1,'LineStyle','--');
    title(forecast.obs{j})
end

set(gcf,'PaperPosition',[0.5 0.5 7.5 10]);
hand=suptitle(['Forecast Observables: Forecasts starting ',forecast.forcsampcell{2}]);
set(hand,'FontSize',16,'FontName','Helvetica');
hand = subtitle(extrsubf(location.savePath,cucd));
set(hand,'FontSize',10,'FontName','Helvetica');
print(gcf,'-dpdf',fullfile(location.savePath,'Forecast Observables'),'-append');

%% Forecast States
AllstateForc=[KFStru.smoothSt(end,:);stateForc]; 
AllstateForc_Small=AllstateForc(:,forecast.location);
numForc=numel(forecast.states);
figure;

for j=1:numForc;
    subplot(ceil(numForc/forecast.plot.cols),forecast.plot.cols,j)
    plot(forecast.forcVec,AllstateForc_Small(:,j),'LineWidth',2)

    axis tight;
%     xlim([forecast.forcVec(1) forecast.forcVec(end)+0.05]);
    title(forecast.names{j})
end

set(gcf,'PaperPosition',[0.5 0.5 7.5 10]);
hand=suptitle(['Forecast States: Forecasts starting ',forecast.forcsampcell{2}]);
set(hand,'FontSize',16,'FontName','Helvetica');
hand = subtitle(extrsubf(location.savePath,cucd));
set(hand,'FontSize',10,'FontName','Helvetica');
print(gcf,'-dpdf',fullfile(location.savePath,'Forecast States'),'-append');